Lindhard and RPA susceptibility computations in extended momentum space in 

electron doped cuprates 
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We present an approximation for efficient calculation of the Lindhard susceptibility x L (c[, ^) in a 
periodic system through the use of simple products of real space functions and the fast Fourier trans- 
form (FFT) . The method is illustrated by providing \ L (cii ui) results for the electron doped cuprate 
Nd2-iCea;Cu04 extended over several Brillouin zones. These results are relevant for interpreting 
inelastic X-ray scattering spectra from cuprates. 
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1. INTRODUCTION 

The dynamic structure factor S*(q, u>) is a useful func- 
tion of momentum and energy introduced by Leon Van 
Hovei, which contains information about density-density 
correlations and their time evolution. Experimentally, 
5(q, uj) can be accessed most directly by Inelastic X-ray 
Scattering (IXS), which has acquired greater importance 
with the advent of powerful synchrotron sources 2 . How- 
ever, since X-rays are strongly absorbed in materials with 
high density, IXS may be suitable mostly for low-Z sys- 
tems. Nevertheless, in the case of heavier elements, re- 
cent studies have shown that if the photon incident en- 
ergy is near an X-ray absorption edge the cross section 
can be enhanced, and the resulting resonant inelastic X- 
ray scattering (RIXS) offers a new window for probing 
both empty and filled electronic states^—. Recent efforts 
to develop a first-principles formulation of the RIXS spec- 
trum explore an interesting hypothesis 6-8 that the RIXS 
cross-section is directly related to S(q, uj), complicating 
effects of the core hole notwithstanding. However, this 
claim remains controversial^ and must be checked by test- 
ing the theory against accurate experimental results^ -. 
Cu-K-edge RIXS for cuprates^iii - — probes the spectrum 
throughout momentum space encompassing many Bril- 
louin zones. Therefore an important theoretical task is 
to produce realistic calculations of the dynamic struc- 
ture factor within the framework of either Many-Body 
Perturbation Theory (MBPT) or Time-Dependent Den- 
sity Functional Theory (TDDFT)— starting from a Lind- 
hard susceptibility representing the response of an un- 
perturbed Kohn-Sham system. In particular, local field 
effects^ - — are known to modify the spectral weight of 
both collective and single-particle excitations in the dy- 
namic structure factor of solids. 

In this study we focus on an approximation to effi- 



ciently calculate from one particle spectral functions the 
Lindhard susceptibility x L (q, w). This approximation 
has successfully described the susceptibility of heavy rare 
earth elements^ and can also reliably describe the X- 
ray inelastic scattering momentum dependency in higher 
Brillouin zones for an energy transfer w where the single- 
particle excitations dominate. As an example we con- 
sider paramagnetic Nd2- 2; Ce a ;Cu04 (NCCO), which has 
a relatively simple, nearly two-dimensional metallic Cu- 
O band near the Fermi level. We identify important fea- 
tures throughout energy-momentum space and we delin- 
eate the specific manner in which Imx (q, uj) decays as 
a function of q. These results enable an assessment of 
the extent to which 5(q, uj) reproduces the RIXS cross- 
section in a cuprate via direct comparison of the theory 
with corresponding experiments in extended regions of 
the momentum space. 

An outline of this paper is as follows. In Section 2, we 
present the relevant formalism. The details of the elec- 
tronic structure methods and the numerical schemes are 
given in Section 3. The theoretical results for Imx L (q, uj) 
are presented and discussed in Section 4, and the conclu- 
sions arc summarized in Section 5. 



FORMALISM 
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In a periodic solid, the susceptibility becomes- 
tensor in the reciprocal lattice vector space G. The 
fluctuation-dissipation theorem relates the dynamical 
structure factor S(q, uj) to the susceptibility via 

s(q, w ) = -ly Im[XG < G(k '" )] 



k,G 



I _ e -huj/kT 



U(q-k-G). (1) 



Thus, IXS experiments do not probe all matrix elements 
of the response xg,c (k, uj), but only the diagonal ele- 



2 



ments XG,G(k, w}2S. If we approximate the susceptibility 
by the bare susceptibility Xg G'(kj w ) thenj^ " 18 ' 21 

Im[ XG , G (k, W )] = -^|M^| 2 f de A v {e)A^ + u). 

(2) 

The matrix elements Mq^ can be expressed in the Dyson 
orbital basis set g v a o 22 ' 23 



--<g v \e^+ G >\ 9fl > 



(3) 



The spectral functions associated with the Dyson orbitals 
are 



A(p,u)^J2\9-(P)\ 2 M^ 



(4) 



with the spectral function A(p, e) expressed in terms of 
the Bloch wave functions instead of plane waves. The ap- 
proximation of Eq. (9) becomes exact when q is large (see 
e.g. Ref. [27]). As already noted above, the asymptotic 
decay of the imaginary part of Lindhard susceptibility as 
a function of q is well described within the present frame- 
work. When q is small, the most significant features of 
the susceptibility are produced by band structure effects, 
which are fully included in our approach. Notably, the 
origin of major peaks in the imaginary part of the sus- 
ceptibility lies in FS nesting . 28 ! 29 Therefore, we expect 
our scheme to produce a reasonable approximation to 
the dynamical structure factor in materials. 



COMPUTATIONAL METHODS 



and 



(5) 



where e v is the excitation energy associated with the 
Dyson orbital g v and 7 is infinitesimally small (see also 
Appendix A). The Dyson orbitals can often be approxi- 
mated reasonably by Bloch orbitals as 

0k,n(r) = exp(ik • r) C G '" exp(-zG • r), (6) 

G 

with the momentum density given by 

|.9k,n(p)| 2 =$>(p-k+G) \C\ 



'k,n|2 



(7) 



G 



The label v = (k, n) is a composite index that codes the 
Bloch wave vector k and the energy band index n. The 
Fourier coefficients C G '™ of the periodic part of the Bloch 
function are labeled by the reciprocal vectors G. In this 
case, the dynamical structure factor at T = becomes 



5(q,w) 



E 



\ya> I G'+G°G» l°G"+G/ 



.k.k'.G'.C'.G 



(8) 

X(5(q + k — k' — G) / deA n (k, e)A m (k', e + lo). 

J —U) 

The dominant part of S*(q, oS) is given by the partial sum 
of the real positive terms G' = G". Next, following 
Wen^ 4 -, we neglect the remaining complex terms because 
the randomness of their phases produces destructive in- 
terferences. A few straightforward algebraic simplifica- 
tions then yield 3 -! 2 ^ 



imx L (q,w) 



de 
2^ 



dp 

— — 3 A(p, e)A(p + q, e + w). 



Thus, our approximation scheme leads to an expression 
for Imx L (q, uj) similar to the free fermion form 2 ^ but 



The Dyson orbitals g v needed for the calculation of 
the spectral function A(p,e), as already noted, can be 
reasonably replaced by the Kohn Sham orbitals obtained 
within the Density Functional Theory (DFT)^ 3 -. For this 
purpose, the DFT band structure calculations in NCCO 
were performed within the Local Density Approximation 
(LDA) using an all-electron, fully charge self-consistent 
semi-relativistic (KKR) method 3 -. The crystal structure 
used for NCCO was body centered tetragonal (space- 
group I4/mmm) with lattice parameters given by Mas- 
sidda et alS^. A self-consistent solution was obtained 
for x = with a convergence of the crystal potential to 
about 1CT 4 Ry. 

To demonstrate our approach in a relatively simple 
but interesting case, we restrict the calculation to a 
single band, namely the copper-oxygen band near the 
Fermi level in NCCO. In particular, the possible con- 
tribution of the Nd f-electrons is neglected by remov- 
ing the f orbital from the basis set after the Nd self- 
consistent potential has been obtained^ 2 - The electronic 
structure shown in Fig. 1(a) has been produced with 
the minority spin part of the self-consistent ferromag- 
netic potential. The doping effects were treated within 
a rigid band model by shifting the Fermi energy to 
accommodate the proper number x of electrons 3 ^—. 
In the electron momentum density (EMD) calcula- 
tions (see Ref. [37] for details), the momentum mesh 
was given by a momentum step (Ap x , Ap y , Ap z ) — 
(l/128a, l/128a,l/2c)/(27r). 3f H° The total number of 
momentum points is 1.54 x 10 8 p within a sphere of radius 
17.6 a.u. 

We show in Fig. [3(a) the calculated band structure 
of NCCO near the Fermi level. The band closest to 
the Fermi level is shown by the red dotted curve and 
is well isolated from other bands. This band ranges from 
— 1.4 eV to 1.9 eV and the integral of the spectral func- 
tion in this energy interval, evaluated with an energy 
resolution of 40 meV, is shown in Fig. [IJb). The two- 
dimensional spectral function A(jp, uS) is calculated by 
neglecting the small k z dispersion in the three dimen- 
sional electronic band structure 41 . Similar EMD results 



for NCCO have been obtained within the LMTO 4 ^. 
The resulting momentum density has the same symme- 
try as the copper-oxygen d x 2_ y 2 —p x , y states in real space 
which form this energy band since the wave function in 
momentum space is the Fourier transform of the wave 
function in real space. Fig. QJb) shows that the low in- 
tensity along the x — y diagonal direction in the 2D-EMD 
map is a signature of d x i_ y i symmetry. Moreover, since 
the radial momentum dependence of an atomic state of 
angular momentum t behaves as p l at small momenta 44 , 
the 2D-EMD intensity at low momenta is from the O- 
2p orbitals, while the Cu-3d x 2_ y 2 orbitals contribute at 
higher momenta! 45 ' 46 This implies that the signal coming 
from the 0-2p states is more visible in the first Brillouin 
Zone while the Cu-3g? states are better seen in higher 
Brillouin zones. We can see from Fig. [IJb) that the 2D- 
EMD intensity is strongly modulated by wave function 
effects, which suggests that the behavior of Imx L (q, w) in 
NCCO in different zones will also be modified by these 
effects; however, our approximation in Eq. (9) neglects 
some interference effects produced by the phases of the 
Fourier coefficients of the Bloch wave functions. 

Equation [S] shows that the Imx L (q, u>) at zero tem- 
perature can be written as a convolution of two spectral 
functions. This Imx L (q, w) captures electron- hole exci- 
tations described by Dyson orbital a 22 ' 23 but does not in- 
clude collective excitations such as plasmons or phonons. 
For efficient Imx L (q, w) calculations, we replace the mo- 
mentum space convolution of A(p, lj) by a simple product 
of spectral functions B(r, lj) in real space given by 

B{r,u) = J A^(p, w )exp( ip r ) . (i ) 

This enables us to take advantage of the fast 
Fourier transform (FFT) efficiency using the convolution 
theorem 4 ^. The advantage of our FFT based method 
can be seen by comparing the computation time of the 
FFT method with the time needed to directly com- 
pute Imx L (q, lj) via Eq. [9] using two matrices of size 
2049 x 2049. The CPU time for the FFT method is 12 
seconds, while the direct computation takes 24 minutes 
on the same machine 4 ^. 



4. RESULTS 

We discuss our results with reference to Figures 2-5. 
In Fig. [5^,, we show — Imx L (q, lj) along high symmetry 
lines as a function of lj. The black part of this figure 
marks the region of zero intensity where no electron hole 
transitions are available. Strong intensity seen near 1 eV 
around (jr, 0) is due to a sort of a Van Hove singularity 
in rmx^q, lj), which is associated with the high energy 
kink or the waterfall effect in the electronic spectrum 4 ^. 
When we compare our Fig. [2^, to the experimental RIXS 
spectrum of overdoped NCCO presented in Ref. we 
find that the experiment is well described by the k re- 
solved joint density of states despite the complicating 
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FIG. 1. (Color online) (a) Band structure of NCCO near 
the Fermi level. The Cu02 band is shown by the red dotted 
line, (b) Calculated integrated spectral function A(p) for an 
isolated Cu02 layer in NCCO. The zone boundaries (for a 
simple tetragonal approximation) are marked by black lines. 
Whites denote large values of A(p), blues small values. The 
A(p) shown here contain contributions only from the Cu02 
band. 




FIG. 2. (Color online) (a) Calculated -Im X L (<l,uj) of NCCO 
along high symmetry lines as a function of transition energy. 
Whites denote the largest -Imx L (q,u), blacks the smallest, 
(b) The integrated value of — Imx (q, lj) over q vs transition 
energy lj. 



effects of the core-hole^. In particular, the features in 
the lowest experimental RIXS band within the energy 
range of lj = 0.5 to lj — 2eV are well reproduced by our 
calculations. The integrated value of — Imx L (q, lj) over 
q, plotted in Fig.[2p, yields the total number of electron- 
hole transitions at a given energy. Since the highest peak 
in Fig. [2b is located at 1.04 eV, we focus on analyzing 
— Imx L (q, lj) at this particular energy in the remainder of 
this article. — Imx L (q, lj) is shown in Fig. [3]for lj = 1.04 
eV over several Brillouin zones marked by yellow lines. 
The first Brillouin zone, located at the center of the fig- 
ure, has the highest intensity. The intensity is seen to 
decrease slowly as q increases, and interesting patterns 
due to d electron wavefunction effects appear in higher 
zones. In the first zone, Fig. [31 some strong peaks are 
present surrounding the zero-intensity hole centered at 
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FIG. 3. (Color online) Calculated -lm X L ((l,uj) of NCCO 
at uj — 1.04 eV plotted over several Brillouin zones. Black 
circles near the F-points indicate regions of zero intensity, as 
in Fig. 2. Index I labels the first zone, indices II, III mark 
zones along the (tt, 0) direction, and IV and V the zones along 
the (tt, tt) direction. 
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FIG. 4. (Color online) (a-d) Contour maps of — Imx L (q, oj = 
1.04 eV) of NCCO in four different Brillouin zones. The sym- 
metry of the maps in zone II and zone III is very similar. 



r with a relatively low intensity appearing at the zone 
corners M. 

Further details of — Imx L (q, uj) are shown in Figs.[4|a)- 
(d), which are blow ups of the four Brillouin zones marked 
by II, III, IV, and V in Fig. [31 The Brillouin zones dis- 
played in Figs. BJa)-(d) show a similar overall pattern 
but modulated with subtle matrix element effects. For in- 
stance, regions of strong intensity spread towards (—Tt, 0) 
in zone II (Fig. @^), but towards (— tt, ±7r) in zone III 
(Fig.0p). The intense (bright) peaks point along one di- 




1 2 3 4 5 6 
q (a.u.) 

FIG. 5. (Color online) Cuts of -Imx L (q,o; = 1.04 eV) taken 
along (n, 0) in Fig. 3 (extended to higher BZs). The high- 
est intensity has been normalized to one. The zone-to-zone 
change in — Imx L (q, uj — 1.04 eV) is mainly an overall de- 
crease in the intensity with smaller changes due to matrix 
element effects. The envelope of — Imx L (q, uj) is fit by us- 
ing a simple Lorentzian shape given by the red dashed line. 
Labels I, II and III correspond to the zone indices in Fig. 3. 
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q (a.u.) 

FIG. 6. (Color online) Cuts of -Imx flPA (q, w = 1.04 eV) (in 
red dashed line) and — Imx L (q, uj) (in blue solid line) taken 
along (7r, 0) in Fig. 3 (extended to higher BZs). The highest 
intensity of each curve has been scaled by taking the highest 
intensity of —Imx L {<l,u) as unity. Labels I, II and III cor- 
respond to the zone indices in Fig. 3. No plasmon peak is 
present in this energy slice uj = 1.04 eV. 



agonal direction in zone IV in Fig. Hfc, but are rotated by 
90 degrees in zone V in Fig. |4ji. 

Figure [5] presents a cut through — Imx L (q, uj) along 
the [100] direction in order to illustrate the decay of 
— Imx L (q, uj) as a function of momentum transfer q. The 
highest intensity has been normalized to unity for ease 
of comparison. Surprisingly, at momenta as large as 6 
a.u. one can still see features with amplitude exceeding 
10 % of the highest intensity (located in the first Bril- 
louin zone). This effect can be explained by the fact 
that d electron particle-hole transitions can involve par- 
ticularly high momentum transfers. We can fit the enve- 
lope of — Imx L (q, uj) by using a simple Lorentzian shape 
. , A v > with fen = 1-89 a.u. 51 Since RIXS has often been 

thought to be related to S(q, uj), it is an interesting ques- 
tion whether a similar decay factor ko is found in RIXS 
experiments. Our results thus provide a new way to test 
the hypothesis that the RIXS cross-section is directly re- 
lated to S*(q, uj). 

We have also obtained the real part of Linhard suscep- 
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tibility \ L by applying Kramers-Kronig relation. How- 
ever, the Linhard susceptibility gives only the response 
of the independent electrons to the external potential. 
In order to estimate the effect of screening effects, one 
can consider the susceptibility within the random phase 
approximation (RPA) given by ^ RPA = x L /(l — K?X L ), 
where V is the Coulomb interaction decaying as ~ 1 /q 2 . 
In this approach, the sharp singularities of x RPA due to 
the denominator give the plasmon modes. Figure 6 il- 
lustrates the corrections are more important when the 
external perturbation is of very long wavelength (i.e. 
q small) 53 . Interestingly when q is of the order of fco, 
— Imx RPA (q, u>) recovers back to the — Imx L (q, w). 



5. CONCLUSIONS 

We have presented a formalism for a first principles 
computation of the Lindhard susceptibility ^(li^) m 
extended momentum space. We have demonstrated a 
tremendous improvement in performance by calculating 
Imx L (q, w) through an approximation involving prod- 
ucts of real space spectral functions B(r, w) and FFTs in- 
stead of using the standard approach involving costly ma- 
trix multiplications. Our theoretical Imx L (q, uj) results 
for the doped cuprate NCCO will allow a detailed com- 
parison with the RIXS experiments, and hence an assess- 
ment of the extent to which Imx L (q, u>) represents a good 
approximation to the RIXS cross section. The present 
work also provides a realistic linear response based start- 
ing point for developing a many-body perturbation the- 
ory of particle-hole excitations within the DFT frame- 



work. 

We are grateful to J. Lorenzana for discussions. This 
work is supported by the US Department of Energy, 
Office of Science, Basic Energy Sciences contracts DE- 
FG02-07ER46352 and lde-sc000709Tl fCMCSN). and ben- 
efited from the allocation of supercomputer time at 
NERSC and Northeastern University's Advanced Scien- 
tific Computation Center (ASCC). It was also sponsored 
by the Stichting Nationale Computer Faciliteiten (NCF) 
for the use of supercomputer facilities, with financial sup- 
port from NWO (Netherlands Organization for Scientific 
Research) . 



Appendix A: Relation between susceptibility and 
spectral function 

We introduce the susceptibility matrix element 



tuu + e„ - e M + 27 ' 



(AI) 



where /(e) is the Fermi function. The term Im[F </ ''j can 
be also written in terms of the spectral function A v . By 
using 



A„(w) = Im[- 



1 



(A2) 



we obtain 



Huj — t v + 17 ' 
Imp"-"] = - f de A„(e)A^e + w) . (A3) 

J —10 
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